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Abstract 

The present paper shows an efficient numerical 
procedure for solving a set of nonlinear partial 
differential equations, specifically the steady 
Euler equations. Solutions of the equations were 
obtained by Newton's linearization procedure, com- 
monly used to solve the roots of nonlinear alge- 
braic equations. In application of the same 
procedure for solving a set of differential equa- 
tions we give a theorem showing that a quadratic 
convergence rate can be achieved. While the domain 
of quadratic convergence depends on the problems 
studied and Is unknown a priori, we show that 
first- and second-order derivatives of flux vectors 
determine whether the condition for quadratic con- 
vergence Is satisfied. The first derivatives enter 
as an Implicit operator for yielding new Iterates 
and the second derivatives Indicates smoothness of 
the flows considered. Consequently flows Involving 
shocks are expected to require larger number of 
Iterations . 

First-order upwind discretization In conjunc- 
tion with the Steger-Warmlng flux-vector splitting 
Is employed on the Implicit operator and a diagonal 
dominant matrix Is resulted. However the explicit 
operator Is represented by first- and second-order 
upwind dlfferenclngs, using both Steger-Warmlng's 
and van Leer's splittings. We discuss treatment 
of boundary conditions and solution procedures for 
solving the resulting block matrix system. With a 
set of test problems for one- and two-dimensional 
flows, we show detailed study as to the efficiency, 
accuracy, and convergence of the present method. 

Introduction 

During the past decades, much progress has 
been made In the development of efficient numerical 
methods for solving unsteady fluid dynamics equa- 
tions to yield time-asymptotic steady solutions. 

The most noted algorithms have been the explicit 
schemes such as the predictor-corrector scheme by 
MacCormackl and more recently the application of 
Runge-Kutta method by Jameson et al., 2 and the 
Implicit ADI schemes such as those due to Briley 
and McDonald 3 and Beam and Warming. 4 Various 
strategies have been utilized to accelerate compu- 
tational efficiencies, e.g., multigrid, 5 vectorl- 
zatlon, etc. 

The present paper proposes a different 
approach wherein the steady-state solutions are 
obtained by solving the steady equations, rather 
than through the time marching approach as 
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mentioned above. We are motivated In this research 
by asking ourselves the following questions: 

(1) Is time marching technique more efficient? 

(2) Can a mathematical basis (even approxi- 
mated) be established to Indicate fast convergence 
of a steady-equations solver? 

(3) How can It be done and what are the dif- 
ferences between the unsteady- and steady-equations 
approaches? 

Obviously the first question Is general and 
the answer Is dependent upon circumstances con- 
sidered. In fact Is Is well known that a space- 
marching approach Is much more efficient than the 
time-marching approach for solving steady, fully 
supersonic problems. Therefore It Is Important to 
Investigate the detailed differences In concept 
and mathematics between the two approaches. We 
note that this research was motivated by the paper 
by MacCormack 6 In which the Idea of using Newton's 
method was mentioned while again unsteady equa- 
tions were used In applications. In this paper we 
attempt to shed some light on the questions (2) and 
(3) together with some applications to various 
problems . 

We apply the Newton method, which Is commonly 
used for Iteratively solving roots of nonlinear 
algebraic equations, to find solutions to the 
steady Euler equations. Consequently quadratic* 
convergence can be achieved provided the Iterative 
sequence Is sufficiently close to the true solu- 
tions. Upwind differencing Is used to approximate 
the derivatives and a diagonal dominant matrix 
system can be obtained, thereby stabilizing the 
Iteration process. Several Investigations® -9 
utilizing Newton's Iteration between two time 
levels and upwind-differencing, have shown great 
success In efficiently arriving time-asymptotic 
steady solutions. It Is noted that the present 
Newton linearization form Is recovered If we let 
At = » In the Implicit unsteady procedures, 
e.g., Beam-Warming's unfactored form. 4 While 
some analytical similarities may exist between the 
steady and unsteady approaches, we believe that 
there are differences In detail which are mani- 
fested In the convergence history. The Iteration 
In general consists of two stages, the Initial 
stage beginning with an arbitrary guess sets out 
to home In the domain of Newton's quadratic con- 
vergence and the second stage locks In this range 
with rapid convergence. 

We show In section 1 a theorem Indicating 
quadratic convergence for the case of differential 
equations. A condition for the domain of quadratic 
convergence, «2» Indicates that whether an 


*Quadrat1c convergence means the error lletl Is 
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approximation lies In G? depends on rate of 
change and smoothness of flow vectors, hence Is 
problem dependent. The choice of spatial differ- 
encing Is of particular Importance for the present 
method and Is given In section 2. Section 3 gives 
treatment of boundary conditions. The algebraic 
system of equations resulting from sections 2 and 
3 Is then summarized; the solution strategies are 
discussed In section 4. Finally we show In sec- 
tion 5 the convergence of calculated solutions by 
comparing with exact solutions In one- and two- 
dimensional problems. The flux splittings by 
Steger and Warming^ and van Leer ^ 2 are employed. 
Pertinent Issues such as accuracy and convergence 
rate are discussed. 

1 . Newton’s Method 

We begin by considering the steady Euler 
equations In conservation-law form and Cartesian 
coordinates : 


¥(U) = D F + D G + OH 

A J ^ 
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(la) 


where ¥ * ,* 2 ¥ 5 } represents five (5) 

partial differential equations, each Is a function 
of the conservative variables U = {p(l ,u,v,w,E)}. 
Here p Is the density, (u,v,w) are Cartesian 
components of the velocity, and E Is the sum of 
specific Internal and kinetic energy. The flux 
vectors are shown below for completeness. 
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In Eq. (la) we choose the notation 
D 


= i- D =i- 

'x ~ 3x ’ y ~ 3y ’ 


•.■fe 


(1C) 


Also perfect gas Is assumed so that flux vectors 
are homogeneous function of U of degree one. 

Assuming U* Is a root (solution) of 
¥(U) - 0, l.e., 


where 6 U » U n+ ^ - U n . The differential operator 
D u * working on $U Is given by 

o/ 5 J* ( U "> = < D X A + °y B + ° 2 C > n (4) 

here the Jacobian matrices are A * D U F, B = D U G, 
and C = D u H.4 

Equation (3) Is the basic equation for deriv- 
ing the new iterate U n +1 and Is, In fact, Iden- 
tical to that derived in Ref. 4 when At * » Is 
taken there. 

For the quasi one-dimensional Euler equations, 

(5) 
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where 
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The Newton procedure yields 


D u ¥ n 4U n = D x (A n • 4U n ) + b" • iu” = -( D X F + S)" 


( 6 ) 


Here we define B as 
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We turn now to study the convergence property 
of the Newton method. Assuming that ¥ Is 
twice differentiable and that the root (solution) 
U* that we are seeking Is a simple root (multi- 
plicity =1). Then we have D U ¥(U*) * 0 and 
D U ¥(U) * 0 for all U In a certain neighborhood 
of U* (see Ref. 10) . 


have 


Using (,) to denote partial derivatives, we 
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*(U*) = 0 


( 2 ) 


Since the system of equations In Eq. (la) Is non- 
linear and coupled, thus we are obliged to use 
approximate method to Iteratively arrive at 
approximate solution to U*. The Newton Iterative 
method having quadratic convergence for solving 
algebraic equations also can be extended to solve 
differential equations. The Newton procedure 
starts by keeping the first-order term In Tavlor’s 
series expansion of ¥(U n+1 ) and forces U n+ ' to 
satisfy Eq. (la) . We have 


D u >? n • «U = - *(0") 


n = 0,1,2,. 


(3) 
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Now expand ¥(tJ ) about *(U ) In Taylor's series 

iyu‘> 

• iee ■ >j) ■ (v >;) <» 


2 


for 1=1,2, P, where V lies In the 

Interval containing U* and U n . Since 

We have, after applying Eqs. (3) and (7), 

‘ " ' “j)(“k ' U k) 

Let c n be the error at the nth Iteration, 
e" E u" - U* 
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Now for maximum vector norms and corresponding 
matrix norms we obtain 

( 12 ) 

where m Is chosen such that 


and denote A = {a^j}, B = {b^j}, C = {c^j}. Sub- 
stituting Eq. (4) In Eq. (8) yields 
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Integrating over an arbitrary volume Ac* = Ax Ay az 
( here ax, Ay, Az are not taken to be grid spac- 
Ings) and making use of the Gauss theorem, we find 

j +1 (a^j dy dz + b^ dx dz + c^ dx dy) 




for all 1 ^ k ^ p 


(13a) 


where the square matrix Q E (q } Is assumed 
Invertible and p Is the number^of equations In C. 
For Euler equations, Eq. (13a) gives, for all k 
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A similar condition for the one-dimensional 
equations Is readily obtained. 
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where (,k) as before denotes differentiation with 
respect to U|<. 

By the mean value theorem the equation becomes 
V/V*) n+l , , n+1 . . c V1 (f) n+1. .A 

e i <*> ♦ (*>) + -isr c j (c) j 
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1, 2, 3. 

Suppose that ¥(U) =0 has a root U* 

In the sphere ^2 = {U:m||u - U*||o> < 1). 

¥ has first and second derivatives such 
(12) Is satisfied. Then for all 


U^S^* n = 1 ,2, . . . 

11m U n = U* , 
n-*» 

the convergence Is quadratic 
* 

U Is unique In 


where 5, n, £, 5 *, V. C « A a. We understood 
that a^j, b^j, and c^j on LHS were evaluated at 
un and a^k, b^j^, and on RHS at V. 

Assuming A<?’ small and 

CjU) * Cj(n) 55 c j(£) etc - 


Proof. Since 


lime 


^ 1, then by Eq. (12) 


lime 




^ 1 
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Hence By Induction It follows for 

n s 1» 2, 3, . ; .. 

Next, repeated use of Eq. (12) gives 




(H) 


Then result (b) follows Immediately as n -» ®. 

Since U°, U 1 , U 2 , . . Is a sequence which 
converges to U*, l.e., U n -» U*, V -> U* as n *♦ », 
there exists a constant C(U*) * 0 such that 
Eq. (12) yields 


After solving a block matrix system (given In 
section 4), we update the approximation by the 
equation 


„n+l • ■ n elt n 

U = U + o)6 U 


(15) 


where the relaxation factor u In general can be 
a constant matrix varying with Iteration. The 
precise variation of o) for which the Iteration 
will remain stable Is not known. A series of 
numerical experiments for various circumstances 
have been performed and their results will be pre- 
sented later. 


Il« n+1 


oo 



n -> » 


The sequence Is, by definition, said to be quad- 
ratlcally convergent. Next we prove the uniqueness 
for a solution In 

Suppose there exists another root V* * U* 

In Q 2 , then *(V*) - For 1 = 1, 2, . . ., p, 

we write 

* * * */* * \ 

U 1 - V 1 = U 1 4. ^(U ) - ^ + ^(V )J 

= U 1 - V 1 + f - v j) 

* ★ 

for V In an Interval containing (U , V ). Hence 
i * *i l * *1 ii * *ii 

h - v ih h - v il + iu - v L^hj< v )| 

Since this Inequality Is true for all 1, 
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a contradiction, hence (d) Is proved. 


We show now the difference between the 
time-marching and the present approaches. Let 
6U n = (i)6U n . Since w Is constant, combination of 
Eqs. (3) and (15) gives 

D • iO n = -(o'? 0 

U 

(16) 


Hence we see that the present Iteration process 
bears no direct relation to the time-marching 
approach since a different Implicit operator hav- 
ing no time derivative term Is Involved. However, 
like the case of time-dependent approach, the 
solution satisfies ¥(U) =0 as 6U -> 0. 

2. Spatial Differencing 

Since there can be Independent choices of 
types of differencing for approximating the Impli- 
cit operator on LHS and explicit operator on RHS 
of the Eqs. (3) or (6) Insofar as stability Is 
maintained, we shall describe them separately. 

The effects of mixed use of differencing on the 
convergence rate and accuracy will be discussed In 
section 5. It will become clear In what follows 
that It Is Imperative to use upwind differencing 
on the Implicit operator. 


From Eq. (14), one can readily estimate the 
number of Iterations required for an Initial error 
to be reduced by 10 _s . 

log 2 109 ^ log (v m | e °lj] 

We shall note that the domain Q 2 I s the 
estimates for the required closeness of approxima- 
tion to the true solution In which Newton's quad- 
ratic convergence holds. It Is function of the 
gradient of the flux vectors with respect to U, 
l.e., Jacoblans, and smoothness of fluxes In the 
problems concerned. 


Now we split the Jacobian matrix A as 

A = A* + A' (17a) 

where A + and A~ are matrices derived from the 
corresponding "+" and split fluxes, F + and 
F” given by Steger and Warming, 11 specifically 



These are the true Jacoblans of relevant split 
fluxes and not Identical to the matrices resulting 
from eigenvalue splitting, as given by: 


In general U Is of course unknown, the 
above conditions, 13(b) or 13(c), must be problem 
dependent and Is difficult to verify beforehand. 

We shall demonstrate through test cases that In 
reality one need not start with close approxima- 
tion. A relaxation procedure Is used to move 
approximation progressively toward the domain Q 2 ; 
thereafter a distinctly faster convergence Is 
followed . 


i 

A" = QA'Q 

where Q Is a similarity matrix diagonalizing A. 4 
The complete expression for A 1 and B* In two- 
dimensional case Is given In Appendix A. 

Now we can write (omitting Iteration Index n 
hereafter) 


4 



D x (A * 4U) = k [A A+4U + A+A 4U] + 0(Ax) (18) 


Here we denote a“ and A + as backward and for- 
ward difference operators. Let 

| A | = A* - A" 

Equation (18) becomes, for the jth grid point, 


D X (ASU) 



•V* A j+i 


-M) 

4U J 

KJ 


(19) 


for ] s 1, 2, . , 3. Substituting of Eq. (19) 

In Eq. (6) gives, for one-dimensional equations, 

D u *6U = T6U (20) 


where T Is a block trldlagonal system {T-| , 

T 2 Tj} and Tj = 1/Ax {-Aj-1. IAjl 

+ AxBj , Aj + i } . Also hereafter 611 denotes {6U 0 , 
6Ui, . . . , 6Uj . . 6Uj, AUj+i); 6U 0 and 
6Uj + i are evaluated at boundaries. 

Following the same procedures In deriving 
Eqs. (12) and (13), we now write their counter- 
parts for the discretized algebraic system, 

It 1 Ila>|d1 scretl zed 4 m (21) 

It Is obvious that whether the condition for quad- 
ratic convergence Is satisfied and fast convergence 
can be realized depends on (1) the property of the 
matrix T resulting from discretization of the 
Implicit operator, and (2) the smoothness of flux 
vector as given In Therefore a well- 

conditioned Implicit operator Is required. Since 
the first-order upwlndlng yields a matrix which Is 
diagonally dominant and hence well-conditioned. 

It also has minimum bandwlth for coupling neigh- 
boring points, hence making boundary conditions 
felt In each Iteration step. Therefore we shall 
use only first-order upwlndlng for the Implicit 
operator while having the liberty of utilizing 
different d 1 f f erenclngs for the explicit operator. 

We note that despite the central differencing 
possesses properties of being simple and higher- 
order accurate with the same bandwlth, It however 
leads to vanishing of diagonal terms and hence 
losing diagonal dominance. Consequently the Iter- 
ative sequence will likely diverge. 

We turn now to the differencing of the 
explicit operator. First-order and second-order 
upwlndlng based on flux splitting are employed on 
uniform grids, we have 

D x F = °X (F+ * F_) (22a) 

and 

°X F± ' k ( A?F± * 2 A ’ A+E?F± ) + 0 [(A*) 1 ” 6 - AX 2 ] 

(22b) 


where E Is the displacement operator, e = 0 and 
1 denote first- and second-order accuracy. Similar 
differencing can be made for D y G and D Z H for 
multidimensional flows. The flux splittings proposed 
by Steger and Warming' 11 and Van Leer 12 are employed 
In our study. They are given here for completeness 
(for two-dimensional case). 

Steger-Warmlng: 

Let F* = 2(y - 1)X* + X* + X* 


F + _ fi— 
h ” 2y 


ufJ + c(X * - xj) 


vf: 


3- F* + CU(Xg - xj) ♦ C 2 (X* + xJ)/2(y - 1) 

(23a) 

where xt = (X, + IX.D/2, 1 = 1, 2, 3, 4 and X. 

= X, = u, X, =u + c, and x. = u - c; c = speed if 
soufid. 3 4 

Now If g| = 2( y - 1)X* + X* + x! and v Is 
substituted for u In x^ above, we have 




uG-j 

VG* + c(X* - xj) 

2 

f- G* + cv(Xg - xj) + c 2 (Xg + \\)/z(y - 1) 

(23b) 


Van Leer: 


F + = F 


F = F 


and for |M| < | 


M > 1 
M < - 1 


(24a) 


r± ' f i 


1 

[(Y - Du ± 2c]/y 
v 

[<Y - Du ± 2c] 2 /2(y 2 - 1) + v 2 /2 


(24b) 


where 

F* = ± pc(M x ± l) 2 /4, M x = U/C . 
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The formula for G~ are obtained by Interchanging 
u and v and elements 2 and 3 in F±. Differ- 
ences of analytical properties of both splitting 
are given In Ref. 12. We shall discuss their dif- 
ferences In calculated results In section 5. 


Now substituting Eqs. (20) and (22) In 
Eq. (6), we have the block trldlagonal system 

T6U = - f (25) 

where f = {fi , f 2 fj fj} and 


AX 


a ’ F j + A 




+ s* 


The system needs to be closed by boundary condi- 
tions, to be discussed In the next section. 

For multidimensional problem, similar upwlnd- 
Ing can be constructed. Assuming equal-spacing 
grids In two dimensions, x = j ax, J = 0, 1 , 2 , . 

. ., 0 + 1 and y = k Ay, k = 0, 1, 2, . . 

K + 1 and a = Ay/Ax. The resulting block matrix 
system Is written as 

M6U = - f (26a) 

where for typewriter ordering, l.e., on each 
constant y, we have 

«U = {«U () ,AU 1 6U r ,. . .,*U K+1 ) 

(26b) 


and each element In 6U and f, 

4U k - {4U 0k’ 4U U 4U Jk 4 Vl,k> • 

k = 0,1,2 K+l 


f = {f lk ,f 2k f jk f Jk } ’ 


k = 1,2,3 K 

(26c) 




2k 


and 


k = 0,1 ,2 K+l 



(28a) 


_A 0k' * A lk* + c * B lk^ A 2k 

~ A lk, * A 2k* + ° l B 2k^ r A 3k 


T k * 5? 


” A J-1 t k |A Jk 1 + a |8 Jk L A J+1 ,k 


" A J-l,k |A 0k ! + * ,B Jk ,f A J+1 ,k 


(28b) 

for k * 1, 2, ...» K. We note that T^ Is a 
block trldlagonal system (each constant y) within 
a large block trldlagonal system M. One can also, 
by altering the typewriter ordering, construct 
block trldlagonal system for each constant-x line. 

There exists many classical Iterative schemes 
for solving a large system; several strategies are 
discussed In Ref. 7. We shall outline possible 
candidates for solving the system, Eq. (26), In 
section 4 after the boundary conditions are 
treated. We note here again conditions on the 
boundary lines (k = 0, K + 1 , j = 0, J + 1 ) must 
be Imposed to render the system solvable. 

3. Boundary Condition 

Boundary conditions are treated Implicitly 
and shown below for various cases. 


The matrix M has the structure 


M = 


Bq T ! B 2 


- B ! T 2 B 3 


B k-1 T k B k+1 


" B K-1 t k B K+1 


(27) 


(a) supersonic Inflow: 

All variables are prescribed and remain 
unchanged at exterior points, 


61 ) = 0 


(29) 


A'F* = o 

(b) supersonic outflow: 

All variables are extrapolated from Interior, 

e.g., 


A U = 0 

a"f + = 0 


(30) 


which again Is of bjock trldlagonal form. The 
element matrices B£ and Tk are found to be 
of block diagonal and block trldlagonal form 
respectively, namely 


(c) subsonic Inflow: 

We prescribe p, u, and v and extrapolate 
p from neighboring Interior points, e.g.. 
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A + p = 0 


(31) 


(d) subsonic outflow: 

The pressure (p) Is fixed and the remaining 
variables are extrapolated from Interior, e.g., 


A’p = 0 
a“u = 0 
A’v = 0 


(32) 


This case Is more Involved and needs some elabora- 
tion. We Illustrate how the boundary terms on LHS 
are combined with the neighboring Interior terms. 
Since p Is fixed, we have 


6pE = 

u 

~ 2 

- 6 p + UA(pU) 


Applying the extrapolation (E 



^*11 a T2 

a l^\ 

a j+i 4U j+i = 

a 21 a 22 

1 

a 23 


1 

y a 31 a 32 

a 33/ 


(° 

a 12 + ua„ 

a 

= 

0 

a 22 + ua 21 

a 



a 32 + ua 31 

a 


(for 1 - D problem) 


f 


\ 


J+l 

13 

23 

33 


- 6p + u6pu 
&pu 


/ 


J+l 


u2a W 


f*A 

u 2 a 21 /2 


<SpU 

l 

u \V 2 J 

J + l 

Wo 


(33) 

Thus we related 6Uj + i with 6Uj and the coeffi- 
cient matrix of 5Uj Is modified by absorbing the 
matrix In Eq. (33). 


chance to get Into the domain $2 given In the 
Theorem shown previously. How fast Is the overall 
convergence rate will then be problem-dependent; 
problem with mixed signs of eigenvalues In general 
will need larger number of Iterations than that 
having only one sign. Furthermore, for a given 
set of Initial guess (U°) and Imposed boundary 
conditions, It Is not possible to estimate a 
priori the domain 1° a nonlinear problem. 

Therefore In practice the entire Iteration sequence 
may consist of two stages. The Initial stage 
starts out an initial guess, which frequently Is 
assumed using some physical judgment, and moves 
slowly to by some suitable relaxation strate- 
gies. We note that one need not use Newton's 
method In this Initial stage and In fact some 
linearly convergent methods may be preferred and a 
better Implicit operator perhaps Is more Important. 
Nevertheless since Newton’s method (Eq. (3)) com- 
bined with relaxation procedure (Eq. (14)) Is 
already losing feature of quadratic convergence, 
hence conceptually It may be thought as some lower- 
order method. Yet this has the simplicity of using 
only one solution procedure and quadratic conver- 
gence (or nearly) is gotten naturally once Q 2 
Is reached. 

For one-dimensional problem, we write the 
block tridiagonal matrix T In terms of lower and 
upper block bidiagonal matrices L and U, l.e., 

LU factorization. 

T = LU (35) 

A description of the procedure Is given In Ref. 13. 
This Is the algorithm that we used to solve one- 
dimensional test problems and may be an optimal 
one considering number of the algebraic operations 
and convergence rate of the overall Iteration 
sequence. However for two-dimensional problems, 

It Is not all that clear which solution scheme for 
a large matrix system, e.g., direct Inversion 
versus Iteration, Is more suitable Insofar as 
efficiency Is concerned. We turn next to outline 
some possible strategies for the system Eq. (26). 

(a) Complete LU Factorization 


(e) solid surface: 


Exactly as above In Eq. (35) we write 


We require the normal component of velocity 
vanishes, l.e., 

V = 0 (34) 

n 

and extrapolate the remaining variables from 
Interior. 

4. Solution Procedure 

With the boundary conditions Implemented, the 
first element of the first row and last element of 
the last row In the block tridiagonal matrix (T 
and M) are essentially eliminated. Since upwind 
scheme provides diagonal dominance, and assuming 
T(M) nonsingular, we can get solution for each 
Iteration. 

Whether this Iteration sequence will converge 
to the true solution (U*) depends on well- 
condltlonedness of the matrix T(M) at each Itera- 
tion, l.e.. If T(M) Is Ill-conditioned, the 
Iteration may diverge quickly and will not have a 


M = LU (36) 

where L and U are of the form 


■ h 




_I 1 r l 


- B I 

L 2 


, u = 

l 2 r 2 



- b k-i 

4 


\ 

1 

l 

u 


The recursive formula for and are given 
In Ref. 13; Inversion and multiplication of block 
matrices are required along each constant-y line. 
The direct Inversion procedure may be too costly 
to perform. In what follows we shall discuss some 
Iterative strategies which can be viewed as m(>l) 
subiterations for solving each Newton’s Iteration 
step. Obviously one need not obtain converged 
solution at each Newton's Iteration which after 
all has not yet produced "true" solution. But the 
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optimal number for m Is not known a priori and 
can be function of many factors, e.g., physical 
problem considered, grid size, etc. In the appli- 
cation of the present method, we chose m = 1 for 
the two-dimensional shock reflection problem, 

(b) Line Relaxation 

If the off-diagonal blocks In Eq. (27) are 
taken to the RHS of Eq. (26a), then we are left 
with a line-decoupled system, l.e., each constant-y 
line Is In effect Independent of each other In so 
far as solving Eq. (26a) Is concerned. Exactly as 
Eq. (23) In form for each k(=l, 2, . . ., K) line, 
we now solve one-dimensional block system, hence 
appropriately called line-relaxation. Since each 
line Is treated as Independent, there Is no 
sequence of sweeping (unlike LU factorization) and 
one can sweep In the same direction for each Iter- 
ation. However there may be advantages alternating 
up and down sweeps In successive subiterations , 
thereby mlmlclng the complete LU factorization. 
Furthermore one may skip odd-numbered lines In one 
sweep and fill In the even-numbered lines In the 
next sweep, l.e., so called the zebra scheme. 
Another possibility Is the Gauss-Seldel updating 
whenever lastest solutions become available, hence 
yielding faster propagation of Information. This 
line Gauss-Seldel relaxation was used In the pres- 
ent paper. 

(c) Point Relaxation 

To further reduce matrix operations, but at 
the expense of slower propagation of Information, 
all off-diagonal blocks In M and T|< are put 
on RHS. An essentially point-decoupled, termed 
point relaxation, system Is resulted. Similar In 
spirit to the zebra scheme, a checkerboard scheme 
can be employed In which points are solved by 
skipping Immediate neighboring points In one sweep 
and filling In those just skipped In the next 
sweep. Again one can also combine Gauss-Seldel 
updating. 

5. Results 

Computational tests of the present method 
were made for one- and two-dimensional flows. 

Four cases of one-dimensional nozzle flows were 
calculated with different Inflow and outflow con- 
ditions. A calculation of oblique shock reflec- 
tion was also made. In all calculations presented 
here uniform grids were used. Inflow conditions 
were also assigned to be the Initial guess for the 
Iterative process. The error Is taken as the sum 
of relative absolute changes of all variables at 
all grid points , l.e., 

error = £ £ III r/ (no. grid points) 

grids 1 l u ll / 

(37) 

so that the convergence of the entire discretized 
system, rather than a single component (variable) 

Is measured. 

In each case calculated we shall present the 
results which were obtained by using combinations 
of (1) first- and second-order dif ferenclngs and 
(2) Steger-Warmlng and Van Leer splittings on the 


RHS explicit operator. Their effects on accuracy 
and convergence will be discussed accordingly. 

(a) One-dimensional nozzle flows 

Two geometries were considered, the area dis- 
tributions are given by: 

Divergent nozzle, 

area = 1.398 + 0.347 tanh (0.8 x -4), 0 < x < 10 

(38a) 

Convergent-divergent nozzle, 

1 1.75 - 0.75 cos (x - 5)*/5 0 < x < 5 

1.25 - 0.25 cos (x - 5 )tt/ 5 5 < x < 10 

(36b) 

We note that since the curvature of the convergent- 
divergent nozzle Is discontinuous at the throat, 
the flow variables (p, p, u) have discontinuities 
In the first derivative. The grid size (Ax) was 
0 . 1 . 

Figures 1 show the calculated pressure dis- 
tributions for fully supersonic flow In the diver- 
gent nozzle, together with exact solution denoted 
by the solid line. Both first- and second-order 
methods gave virtually the same results. It Is 
evident that a truly quadratic convergence was 
obtained using the first-order differencing (error 
was reduced by 10-1* In seven Iterations); here 
g> Is the relaxation factor appearing In Eq. (15). 
However the second-order differencing appeared to 
achieve quadratic convergence Initially, but 
changed to much slower rate. This peculiar behav- 
ior was not understood and Is currently under 
study. 

We now present the case with a shock, l.e., 
the outflow Is subsonic, In Fig. 2. The first 
order upwlndlng produced monotonic behavior; the 
Steger-Warmlng splitting smeared the shock more 
than the Van Leer splitting did, but was about 
twice as fast In convergence rate. The difference 
In convergence rate could be attributed to the 
fact that the Implicit operator remained unchanged 
(l.e., Steger-Warmlng splitting) the explicit 
operators used either splitting. The to was set 
equal to 0.5 Initially and Increased according to 
the equations: 

o> = a> + 0.1 n/n Q If M0D(n,n Q ) = 0 

and 

a> - ml n ( co , 1 .0) 

where n 0 Is the Interval of Iterations In which 
w Is kept fixed. The formula was arbitrarily 
chosen; better convergence can be gotten when the 
relaxation process Is better understood. The 
second-order differencing produced much sharper 
shock, but with unwanted oscillations at both 
upstream and downstream of the shock. The conver- 
gence rate was about half of the first-order 
method. 
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Next we looked at flows In the convergent- 
divergent nozzle. Figures 3 show the fully sub- 
sonic case. Both splitting gave virtually 
Identical results except differences In conver- 
gence; second-order differencing was clearly 
superior In accuracy and also nearly equal In 
convergence rate, contrary to the case In Fig. 1. 

Figures 4 display the case Involving a shock; 
the first-order method clearly lacked spatial 
resolution. The Van Leer splitting gave much bet- 
ter results In the supersonic region but over- 
predicted In the Inflow subsonic region. The 
second-order method obviously showed much Improved 
accuracy, but at the expense of slower convergence 
rate. The dispersive error near the shock was 
observed In both splittings; the discontinuity at 
the sonic point In the Steger-Warmlng splitting 
was evident and did not appear In the Van Leer 
splitting, which was designed to yield smooth 
transition of splitting at sonic point. 


explicit operator, combined with the first- and 
second-order dlf ferenclngs . One-dimensional nozzle 
flows and two-dimensional shock reflection problem 
were calculated, specific emphases were placed on 
the convergence rate and accuracy. 

We are continuing to make In-depth Investiga- 
tion of the present method, specially as to the 
Increase of efficiency and application to more 
complex problems. A high resolution scheme has 
been Implemented 1 * and the results will be pre- 
sented elsewhere. 
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Appendix A 


Finally we turn to the problem of regular 
reflection of a two-dimensional oblique shock wave 
with free stream Mach number 2.9 and shock angle 
29°. The computation domain was divided uniformly 
with 61 by 21 grids over 0.0 < x < 4.0, and 
0.0 < y < 1.0, the Steger-Warmlng flux splitting 
was used for calculating this case. Figure 5 
gives the pressure distributions at the wall, using 
the first-order upwlndlng. The shocks evidently 
were smeared to a large extent. While the static 
pressure showed monotonic property, the total 
pressure (or a measure of entropy) displayed over- 
shoot followed by a recovery. The convergence 
history showed a drastic reduction In error after 
an Initial stage. Figure 6 shows the results 
obtained by the use of second-order upwlndlng 
(boundary conditions remained first-order approxi- 
mation). Here we see an overexpansion In pressure 
just ahead of shock foot, followed by a smooth 
Increase. The spatial resolution was clearly much 
Improved over the first-order method. The over- 
expansion was also accompanied by the Increase of 
total pressure, not seen In the first-order method. 
The convergence history again showed two distinct 
stages, confirming the theorem and speculation 
given In section 1. The convergence history for 
the error to reduce lO- 1 ^ roughly the same 
for both first- and second-order dlf ferenclngs . 

The Mach contours are depicted In Fig. 6(d). 


We show the true Jacoblans which are derived 
from the split flux vectors given by Steger and 
Warming 11 In what follows. 

As defined previously 


F = F + * F", 
A = 3F/3U, 

A 4 = 3F 4 /3U, 

|F| = F + - F". 

Hence [ F | - F for u > 
u < -c. 

Let 

I A | = 

hence 


G = G + + G 
B = 3G/3U 

(AT) 

B” = aG /au 
| G | = G + - G" 
c and | F | = -F for 


A - A 


( A2) 


A 4 = (A ± | A | ) /2 (A3) 

From Eq. (Al) we have 


6 . Concluding Remarks 

We summarize that an efficient method has 
been proposed to find steady solution by solving 
the steady Euler equations. The Iterative scheme 
necessary for solving a coupled nonlinear system 
was based on Newton's method and an updating 
(relaxation) procedure. A theorem was shown that 
the Newton method applied to the system of differ- 
ential equations gave a quadratic convergence pro- 
vided the approximation to the true solution Is 
sufficiently close. The first and second deriva- 
tives of flux vectors were shown to define the 
domain of the quadratic convergence (s^)- Upwind 
differencing was necessary to construct the 
Implicit operator for providing the Iteration 
sequence. We employed the true Jacoblans resulting 
from the Steger-Warmlng splitting In approximating 
the Implicit operator. Both Steger-Warmlng 1 s and 
Van Leer's splittings however were used on the 


|A| - ^u 1 - fu (F+ - F "> ■ {i 13 } (A4) 


Since the fluxes F and G are split according to 
Eqs. (23) and (24), a straightforward substitution 
and differentiation yields, after tedious algebra, 
for | u | < c: 


*11 = “ E * ° 5 X 2T L 


a 12 » ± 13 -ou , 0 = (y - 1 )/y 


*13 " - aV 


a 


14 


= a 
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a ?1 = + Bu 2 + 2au(-e + q 2 /2) , q 2 = u 2 + v 2 


a^ 2 = 2{± Bu + c/y} - 2aU 


3. Briley, W.R. and McDonald, H., "Solution of 
the Multidimensional Compressible Navler-Stokes 
Equations by a Generalized Implicit Method," 
Journal of Computational Physics . Vol . 24, 

No. 4, Aug. 1977, pp. 372-397. 


a 23 = " 2aUV 


a 24 = 2otU 


a 


31 


+ Buv + avl-e 



4. Beam, R.M. and Warming, R.F., "An Implicit 
Factored Scheme for the Compressible Navler- 
Stokes Equations," AIAA Journal . Vol. 16, 

No. 4, Apr. 1978, pp. 393-402. 

5. Jameson, A. and Baker, T.J., "Multigrid Solu- 
tion of the Euler Equations for Aircraft 

( A5) Configurations," AIAA Paper 84-0093, Jan. 1984. 


a 32 = ± Bv - auv 


6. MacCormack, R.W., "Current Status of Numerical 
Solutions of the Navler-Stokes Equations," 

AIAA Paper 85-0032, Jan. 1985. 


a 33 " 1 6U + 


c/y - aV 


a 


34 = 


aV 


a., = E[o(3u 2 + v 2 ) + 3c]/2 

4 I 


7. Chakravarthy, S.R., "Relaxation Methods for 
Unfactored Implicit Upwind Schemes," AIAA 
Paper 84-0165, Jan. 1984. 

8. Thomas, J.L. and Walters, R.W., "Upwind 
Relaxation Algorithms for the Navler-Stokes 
Equations," AIAA Paper 885-1501, July 1985. 


- {q 2 [± 6u + c/y] + 2c[u 2 + c 2 /(y - 1)]} 
a 42 = B[± q 2 /2 + u 2 ] + 3uc/y - u[t»(3u 2 + v 2 ) + 3c]/2 
a 43 = v[± Bu + c/y] - v[o(3u 2 + v 2 ) + 3c]/2 


9. Walters, R.W. and Dwoyer, D.L., "An Efficient 
Iteration Strategy for the Solution of the 
Euler Equations," AIAA Paper 85-1529, July 
1985. 

10. Dahlqulst, G. and Bjorck, A., Numerical Methods 
(translated by N. Anderson), Prentice-Hall, 
1974. 


a.. = o(3u 2 + v 2 )/2 + 3c/2 

44 

where "+" and "-" correspond to u > 0 and u < 0. 
Similar expression for |B| can be obtained by 
Interchanging rows 2 and 3, columns 2 and 3, and 
u and v. 
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NEWTON METHOD-FIRST ORDER UPWIND 



NEWTON METHOD-SECOND ORDER UPWIND 



Figure 3. - Static pressure and convergence history for convergent-divergent nozzle, 
Moo- 0.2, P exjt /P t - 0.9188. 


NEWTON METHOD-FIRST ORDER UPWIND 



NEWTON METHOD-SECOND ORDER UPWIND 



Figure 4. - Static pressure and convergence history for convergent-divergent nozzle, 
Moo- 0.2395; P exjt /P t = 0.84. 
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